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Preface 


A  common  practice  in  reliability  engineering  is  to 
utilize  available  component  sample  data  to  derive  a  point 
estimate  of  the  reliability  of  each  component.  This  thesis 
is  a  continuation  of  previous  work  accomplished  at  the  Air 
Force  Institute  of  Technology  on  reliability  estimation. 

Its  purpose  is  to  indicate  the  manner  in  which  a  modified 
Double  Monte  Carlo  technique  can  be  utilized  to  derive 
confidence  interval  estimations  of  system  reliability  based 
on  sample  component  data. 

I  wish  to  express  my  sincere  appreciation  to 
Professor  Jon  R.  Hobbs  for  his  guidance  and  expertise 
toward  completion  of  this  study  and  to  Professor  Albert  H. 
Moore  for  suggesting  this  topic  and  his  encouragement  and 
guidance  throughout  its  development. 
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Abstract 

A  digital  computer  technique  is  developed,  using  a 
modified  Double  Monte  Carlo  simulation,  which  determines 
lower  confidence  limits  for  system  reliability  based  on  com¬ 
ponent  test  data.  This  test  data  is  assumed  to  have  failure 
times  which  are  distributed  according  to  a  known  two- 
parameter  Weibull  probability  distribution.  The  first  step 
of  the  modified  Double  Monte  Carlo  technique  is  to  randomly 
generate  these  failure  times  using  the  true  shape  and  scale 
parameters  of  each  component.  The  component  distribution 
shape  and  scale  parameters  are  then  estimated  by  the  method 
of  maximum-likelihood  from  these  component  failure  times. 

The  second  step  of  this  technique  is  to  reestimate  these 
component  shape  and  scale  parameters  using  generated  samples 
whose  failures  have  the  same  distribution  and  parameters  as 
the  estimated  ones  and  the  same  number  of  observations  as 
the  original  test  data.  The  method  of  maximum-likelihood 
is  again  used  to  estimate  these  component  parameters.  These 
twice  estimated  parameters  are  then  substituted  into  the 
reliability  equation  to  obtain  the  maximum- 1 ikel ihood  esti¬ 
mator  for  the  component  reliability.  The  estimated  bias  in 
this  estimator  is  subtracted  to  yield  an  approximate] y 
unbiased  estimator  of  component  reliability.  A  given  number 
of  component  reliabilities  are  obtained  and  used  with  the 
median  rank  values  to  construct  an  empirical  distribution 
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for  each  component.  The  desired  number  of  estimates  are 
then  sampled  from  these  distributions  to  obtain  a  system 
reliability.  Using  this  technique,  where  the  distribution 
or  joint  distribution  of  the  estimators  is  unknown,  a 
Monte  Carlo  simulation  is  run  for  four  hypothetical  systems 
consisting  of  as  many  as  five  components.  Since  the  true 
reliability  is  known,  it  can  be  determined  if  the  desired 
confidence  intervals  contain  the  true  system  reliability. 

The  result  is  a  measure  of  the  effectiveness  of  the  modified 
Double  Monte  Carlo  technique. 
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A  MODIFIED  DOUBLE  MONTE  CARLO  TECHNIQUE  TO  APPROXIMATE 
RELIABILITY  CONFIDENCE  LIMITS  OF  SYSTEMS  WITH 
COMPONENTS  CHARACTERIZED  BY  THE 
WE I BULL  DISTRIBUTION 

i 

i 

< 

I.  Introduction  j 

j 

Problem  Statement 

In  the  Air  Force  today,  there  is  a  tendency  towards 
increasing  complexity  of  systems  which  correspondingly  makes 
the  achieving  of  high  orders  of  reliability  more  difficult. 

Because  of  this  increasing  complexity  and  the  concern  over 
the  ever  increasing  Defense  budget,  these  reliabilities  must 
be  evaluated  by  using  less  testing  and  more  efficient  methods 
of  reliability  estimation.  As  a  result,  much  effort  has 
been  made  in  establishing  methods  for  predicting  reliabili¬ 
ties  of  complex  systems  from  their  component  test  data. 

Since  the  exact  reliability  of  a  component  can  not  be 
measured  directly,  it  must  be  estimated.  These  estimates 
can  vary  considerably  in  their  accuracy;  therefore,  it  is 
necessary  that  limits  be  attached  to  their  probable  range. 

When  these  limits  are  determined  to  a  desired  degree  of 
confidence,  reliability  confidence  limits  result.  It  is 
meaningless  to  say  that  a  system  is  90%  reliable,  after 
testing  its  components,  without  some  type  of  confidence 
limit.  The  purpose  of  this  thesis  is  to  determine  these 
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reliability  confidence  limits  based  on  test  data  about 
individual  components  that  comprise  the  system. 

Background 

In  1960,  Orkand  (Ref  14)  used  a  Monte  Carlo  method 
for  determining  lower  confidence  limits  on  system  reliability. 
His  study  showed  the  limitations  in  using  only  point  esti¬ 
mates  of  component  reliability  to  determine  a  point  estimate 
of  the  system  reliability.  Bernhoff  (Ref  3)  explored  the 
problems  of  applying  analytical  approaches  to  establishing 
system  reliability  confidence  limits.  He  concluded  that 
system  confidence  limits  can  be  obtained  analytically  if 
all  components  of  the  system  have  the  same  mathematical  form 
for  reliability.  If  the  system  reliability  is  a  function  of 
two  or  more  dissimilar  mathematical  expressions,  the  only 
practical  way  of  finding  the  approximate  reliability  distri¬ 
bution  is  to  use  a  Monte  Carlo  technique  (Ref  3:49). 

Levy  (Ref  8)  developed  a  Monte  Carlo  technique  where 
components  were  subjected  to  life  tests  and  the  mathematical 
model  for  component  failures  was  assumed  to  be  of  the 
Exponential,  Normal,  Lognormal,  Gamma,  or  Weibull  distribu¬ 
tion  . 

Moore  (Ref  13)  extended  the  Monte  Carlo  methods  to 
include  those  cases  where  the  joint  distribution  of  the 
estimators  is  unknown.  He  was  able  to  estimate  their  joint 
distribution  by  use  of  a  technique  he  developed  called 
Double  Monte  Carlo. 
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Moore  and  Levy  (Ref  9)  designed  a  digital  computer 
technique  for  obtaining  system  reliability  confidence  limits 
where  the  system  component  failures  exhibited  different 
probability  densities.  It  was  assumed  that  the  location 
parameter  is  zero  or  known,  that  the  shape  parameter  is 
known  and  that  the  probability  density  function  is  given. 

With  these  assumptions,  the  exact  distribution  of  the  maxi¬ 
mum  likelihood  estimator  can  be  determined. 

Lutton  (Ref  10)  assumed  that  the  component  life  distri¬ 
butions  were  known.  He  gen'  -ated  reliability  samples  using 
the  original  Levy  (Ref  8)  technique  and  the  asympotic  dis¬ 
tribution  of  parameter  estimates  of  the  Weibull,  Gamma  and 
Logistic  density  functions. 

Lannon  (Ref  7)  established  reliability  confidence 
intervals  for  the  Weibull  distribution  using  a  bivariate 
analysis  with  the  scale  and  shape  parameters  unknown. 

Most  recently,  Snead  (Ref  18)  studied  reliability 
estimates  using  the  Weibull,  Gamma,  and  Logistic  distribution 
with  the  property  that  the  reliability  estimators  are 
asympotically  normal.  This  study  required  large  sample 
sizes  to  perform  a  univariate  analysis  using  just  the  relia¬ 
bility  parameter.  Putz  (Ref  15)  further  developed  this 
univariate  technique  by  simplifying  the  bivariate  analysis 
of  Lannon  (Ref  7),  by  mapping  the  shape  and  scale  parameters 
onto  the  reliability  parameter  and  reducing  the  sample  sizes 
of  the  component  test  data. 
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Rice  (Rex'  16)  incorporated  the  asymptotic  normality 
properties  of  the  binomial  distribution  using  a  Monte  Carlo 
technique  for  estimating  lower  confidence  limits  of  system 
reliability.  He  also  incorporated  a  technique  developed  by 
Gatliffe  (Ref  4),  which  substituted  equivalent  failures  in 
his  component  failure  estimations,  when  the  component 
exhibited  no  failures. 

Darrel  Thoman ,  Lee  Bain,  and  Charles  Antle  (Ref  19) 
showed  that,  assuming  the  shape  and  scale  parameters  are 
unknown,  location  parameter  is  zero,  in  the  two-parameter 
Weibull,  that  the  distribution  of  the  maximum  likelihood 
estimator  of  reliability,  depends  only  on  the  true  system 
reliability  and  the  sample  size. 

Reliability 

The  reliability  of  a  system  is  defined  as  the  proba¬ 
bility  that  the  system  will  be  operating  at  some  specified 
time  (t)  under  specific  conditions.  If  T  is  the  time 
to  failure  or  life  length  of  a  system  or  component,  the 
reliability  at  time  t  or  R(t)  is  given  by  R(t)=P(T>t)  , 

where  P  means  "probability  of".  The  system  reliability  is 
dependent  on  each  component  reliability  and  the  system  con¬ 
figuration.  For  example,  if  the  components  are  connected 
serially,  the  failure  of  any  one  of  them  will  cause  the 
system  to  fail.  If,  on  the  other  hand,  the  components  are 
connected  in  parallel,  the  system  will  fail  only  if  all  com¬ 
ponents  or  a  specified  number  of  them  fail. 


4 


The  Weibull  Distribution 


The  Weibull  probability  density  function  (pdf)  was 
originally  developed  and  used  by  Waloddi  Weibull  in  1939  in 
a  study  of  the  phenomenon  of  rupture  in  solids.  Since  that 
time,  the  Weibull  pdf  has  been  found  to  be  useful  for  appli¬ 
cation  in  lifelength  and  reliability  testing  for  many  mechan¬ 
ical  and  electronic  components.  It  is  frequently  assumed 
that  electronic  components  are  distributed  according  to  the 
Exponential  distribution.  Zelen  and  Dannemiller  (Ref  20:36) 
have  pointed  out  that  the  Exponential  distribution  is  gen¬ 
erally  not  a  robust  approximation  to  the  Weibull,  especially 
if  the  shape  parameter  is  greater  than  1. 

The  Weibull  density  function  is  defined  as 


f  (t  ;k,  0  , c)  =  Ml_c) -  e  V  9  '  k,02O;  c£t  (1) 

0 

=  0  elsewhere 

where  0  is  the  scale  parameter,  k  is  the  shape  parameter, 
c  is  the  location  parameter  and  t  is  the  time.  The  scale 
parameter  affects  the  dispersion  of  the  random  variable  t 
about  its  mean.  The  shape  parameter  determines  whether  the 
hazard  function  is  increasing,  decreasing,  or  time  invariant; 
while  the  location  parameter  determines  the  origin  (or 
guaranteed  life).  The  Weibull  distribution  can  be  used  to 
model  the  exponential  density  function,  if  the  shape  para¬ 
meter  is  1.  It  can  also  approximate  the  Normal  distribution 
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and  the  Raleigh  distribution  when  the  shape  parameter  is 
scaled  at  3.5  or  2.0  respectively.  Figure  1  shows  the 
flexibility  of  the  Weibull  to  model  Exponential,  the  Normal 
and  the  Raleigh  distributions. 

Assumptions 

This  thesis  will  concern  itself  with  systems  which 
have  components  described  by  the  Weibull  pdf.  It  is  assumed 
that  the  components  have  previously  been  determined  to  be 
Weibull  or  that  the  Weibull  pdf  adequately  models  the  com¬ 
ponents  in  the  system. 

It  is  assumed  that  the  components  of  the  system  being 
analyzed  fail  independently.  That  is,  failure  of  a  given 
component  within  a  system  does  not  depend  upon  either 
failure  or  successful  operation  of  any  other  component. 

The  location  parameter  is  assumed  to  be  known  and  will 
be  set  equal  to  zero  for  all  cases. 

It  is  assumed  that  all  components  have  been  life 
tested  and  that  all  unknown  parameters  of  the  life  distri¬ 
bution  have  been  estimated  from  the  data. 

Confidence  Limits 

It  is  known  that  statistical  estimates  are  more 
likely  to  be  close  to  the  true  value  as  the  sample  size 
increases.  Thus,  there  is  a  close  correlation  between  the 
accuracy  of  an  estimate  and  the  size  of  the  sample  from 
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TIME  *10l 

Fig  1.  Weibull  Distribution  with  Shape 
Parameters  of  1,  2,  and  3.5 

which  it  was  obtained.  To  this  extent,  in  order  to  obtain 
a  100  percent  confidence  or  certainty  that  a  measured 
statistical  parameter  coincides  with  the  true  value,  an 
infinitely  large  sample  size  or  infinite  interval  is 
required. 

When  the  estimate  of  a  parameter  is  obtained  from  a 
reasonably  sized  sample,  it  may  be  logically  assumed  that 


the  true  value  of  that  parameter  will  be  somewhere  in  the 
neighborhood  of  the  estimate.  Therefore,  it  is  more  mean¬ 
ingful  to  express  statistical  estimates  in  terms  of  a  range 
or  interval  with  an  associated  probability  or  confidence 
that  the  true  value  lies  within  such  an  interval,  than  to 
express  them  as  point  estimates.  This  is  in  fact  what  is 
accomplished  when  confidence  limits  are  assigned  to  point 
estimates . 

Confidence  intervals  around  a  point  estimate  have  a 
lower  confidence  limit  L  and  an  upper  confidence  limit  U 
For  example,  if  it  is  desired  to  calculate  the  confidence 
limits  for  a  probability  of  90  percent,  this  means  that  in 
approximately  90  out  of  100  cases  the  true  value  will  lie 
within  the  calculated  limits,  and  approximately  10  cases 
will  lie  outside  these  limits.  Distinguished  from  this 
confidence  interval  is  a  confidence  level  where  it  is 
assured  that  at  some  given  level,  say  10  percent,  the  true 
value  lies  within  the  calculated  limits.  If  it  is  desired 
to  increase  the  confidence  level  to  99  percent,  so  that  in 
99  percent  of  the  cases  the  true  value  would  lie  within  the 
confidence  limits,  the  confidence  interval  around  the  point 
estimate  would  become  much  wider--or  a  much  larger  sample 
size  must  be  used  for  the  point  estimate. 

Objective 

Compare  the  relative  accuracy  and  utility  of  the 
modified  Double  Monte  Carlo  technique  for  finding  system 
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reliability  confidence  limits  with  previously  determined 
methods  using  the  standard  Double  Monte  Carlo  method  and 
the  Univariate  and  Bivariate  asymptotic  distributions. 

Approach 

There  must  be  a  known  mathematical  relationship  between 
the  system  reliability  and  component  reliabilities  (Ref  13: 
459).  Hence,  the  assumption  that  the  components  fail  inde¬ 
pendently  of  one  another.  This  allows  use  of  standard 
formulas  to  calculate  system  reliabilities  whether  the  com¬ 
ponents  are  connected  in  series  or  in  parallel.  Orkland 
describes  procedures  to  apply  if  there  is  a  dependence  rela¬ 
tion  among  system  components  (Ref  8:7-8). 

A  modified  Double  Monte  Carlo  technique  will  be  used 
to  determine  reliability  estimates  and  the  associated  confi¬ 
dence  intervals  of  specific  component  networks.  The  unique 
feature  of  this  method  is  that  the  distribution  (joint 
distribution)  of  the  estimator(s)  for  the  parameter(s)  can 
be  unknown  (Ref  13:461).  The  Monte  Carlo  method  uses  random 
sampling  to  investigate  the  solution  of  deterministic  or 
stochastic  problems.  In  essence,  one  inserts  a  set  of  random 
inputs  from  a  specific  probability  distribution  and  solves 
the  problem  for  each  set  of  inputs  to  obtain  a  random 
sample  of  outcomes  (Ref  13:459).  The  steps  to  this  method 
are  as  follows: 
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1.  Calculate  the  maximum  likelihood  estimators  of 
the  component  shape  and  scale  parameters  from  the  component 
failure  times.  Do  this  for  all  components. 

2.  Samples  are  then  generated  whose  failure  have  the 
same  distribution  and  parameters  as  the  estimated  ones  and 
the  same  number  of  observations  as  the  original  test  data. 
These  parameters  are  also  estimated  by  maximum  likelihood. 

3.  Substitute  the  sample  parameter  estimates  into  the 
component's  reliability  equation  to  derive  a  sample  component 
reliability . 

4.  Repeat  steps  1-3  to  form  an  empirical  distribution 
of  component  reliability  estimates  for  each  component. 

5.  Generate  a  point  sample  compone;.t  reliability  for 
each  component  from  (4). 

6.  Calculate  a  point  sample  of  system  reliability 
from  the  point  sample  of  component  reliability. 

7.  Repeat  steps  5-6  to  obtain  the  desired  number  of 
samples  of  the  system  reliability. 

8.  Order  the  system  reliability  samples  to  obtain  a 
sample  cdf  of  system  reliability  from  which  the  confidence 
limits  are  determined  at  a  given  confidence  level. 

The  main  advantage  of  the  Monte  Carlo  method  is  its 
usefulness  in  solving  very  complicated  problems.  The  disad¬ 
vantage  is  in  its  slow  convergence  and  resulting  higher 
computer  costs  (Ref  13:459). 
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I I .  Theoretical  Development 


Estimators 

An  estimator  is  a  rule  which  tells  how  to  calculate 
the  estimate  of  something,  based  upon  information  contained 
in  the  sample  (Ref  2:161).  One  type  of  estimator  is  a  point 
estimate  where  a  single  number,  representing  the  estimate, 
may  be  associated  with  a  point  on  a  line.  An  example  of  a 
point  estimate  is  the  sample  mean  y  ,  where 
n 

£  yi 

y  =  — - -  .  This  estimator  of  the  population  mean,  y  , 

explains  exactly  how  the  actual  numerical  value  of  the  esti¬ 
mate  may  be  obtained  once  the  sample  values  y2»yo.--*»yn 
are  known. 

A  different  problem  with  multiple  parameter  estimation 
may  be  stated  as  follows:  Assume  that  f ( t ,  0 ^  ,  0g , • • • , 6^) 
is  a  density  with  k  unknown  parameters  and  let 
t, ,t0,...,t  be  a  random  sample  of  size  n  .  In  the  case 
of  the  Weibull  pdf,  there  are  three  parameters  k,0,c  , 

with  tlft2,...,tn  being  sample  failure  times.  The  problem 
is  to  find  a  function  of  the  observed  samples  which  may  be 

represented  by  k( tj , t2 , . • • , tn ) ,  0(t1>t2 . tn),  cCt^.tg, 

...,tn)  ,  such  that  the  distribution  of  these  functions 

will  approximate,  as  closely  as  possible,  the  true  values  of 
the  parameters.  Each  of  these  functions  will  be  an  estimator 
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of  the  true  value  and  will  be  denoted  as  k,  0,  c  respec¬ 
tively.  (Note:  In  this  thesis  the  location  parameter  c 
is  always  assumed  to  be  zero. ) 

Method  of  Maximum  Likelihood 

A  very  general  method  of  finding  point  estimates  of 
parameters  is  the  maximum  likelihood  estimator  (MLE).  This 
method  selects  those  values  of  the  parameters  which  maximize 
the  probability  or  the  joint  density  (the  likelihood)  of 
the  observed  sample  (Ref  11:302).  The  likelihood  function 
can  be  defined  as  follows:  Let  y1,y0,...,y  be  sample 
observations  taken  on  corresponding  random  variables, 

Y-^Yg . Yn  .  Then  if  Y1>Y2,...,Y  are  discrete  random 

variables,  the  likelihood  of  the  sample,  L  ,  is  defined 

to  be  the  joint  probability  of  y-^y.^ . Yn  • 

Y1,Y2>...,Yn  are  continuous  random  variables,  the  likeli¬ 
hood,  L  ,  is  defined  to  be  the  joint  density  evaluated  at 

y^.yg . yn  (Ref  11:303).  For  the  Weibull  density,  the 

likelihood  function  can  be  represented  as: 

L  =  f(t1.t2. • • • »tn;  k , 0 , c  =  0) 

Since  there  is  an  assumption  that  the  failure  times  are 

n 

independent  then  L  =  tt  f(t.;  k,8) 

i  =  l 

/S  A 

It  is  desired  to  determine  k,0  ,  the  MLE  for  k 

and  0  respectively,  the  procedure  used  in  this  thesis  is 
as  follows: 


12 


1.  Take  the  partial  derivatives  of  the  natural 
logarithm,  with  respect  to  each  parameter. 

2.  Set  these  derivatives  equal  to  zero. 

3.  Solve  simultaneous  equations  for  the  values  of  the 
parameters . 

Since  the  Weibull  pdf  is 


f(t;  k,e,c=0) 


then 


k, 0 , tao 


T  ,  n.-kn,  n  ,  k-lv  „  n-k  ”  .  k. 

L  =  k  0  (  7r  t.  )exp(-0  E  t .  ) 

i-1  i-1  1 


and 


n 

In  (L)  =  n  In  k  -  nk  In  0  +  (k-1)  I  In  t. 

i=l  1 


(2) 


(3) 


-k  n  k 
6  I  t. 

i=l  1 


(4) 


Taking  the  partial  derivative  of  In  (L)  with  respect  to 
0  and  k  , 


=  “nk/6  +  k9~k  1  ^  t.k  =  0 

30  ,=1  i 


(5) 


3  In  L 
3k 


=  n/k  -  n  In  0  +  E  In  t. 

i=l  1 


.  n  . 

-  0_K  E  t,  In  t.  =  0  (6) 

i=l  1  1 
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Denoting  the  solutions  by  k  and  0  ,  Eqs  (5)  and  (6) 

may  be  rewritten  as: 


0 


n 

Z 

i=l 


V 


(7) 


-  -c  n  C  n 

k  =  n/[l/0K  E  t.  In  t.  -  E  In  t.]  (8) 

i=l  1  1  i=l  1 


A 

Eqs  (7)  and  (8)  are  two  equations  with  two  unknowns,  k 
and  0  .  Simultaneous  solution  of  these  equations  by  an 

iterative  procedure  developed  by  Harter  and  Moore  (Ref  5) 
yields  the  maximum  likelihood  estimators  of  k  and  0 


Component  Reliability 

The  reliability  function,  R(t)  ,  may  be  expressed 

as : 


R(t) 


f (x)dx 


t 


(9) 


where  x  is  dummy  variable  of  integration  and  f(t)  repre¬ 
sents  a  failure  density  function.  However,  if  the  para¬ 
meters  of  f(t)  are  unknown  and  must  be  estimated  from 
data  samples,  the  reliability  function  itself  must  be 
expressed  as  an  estimate;  thus 

R(t)  -  g(t;  0^02 . 6k) 
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Since  each  component  and  their  associated  parameters 
are  characterized  by  the  Weibull  distribution,  the  relia¬ 
bility  function,  R(t)  ,  is  determined  using  Eqs  (2)  and 
(9)  as: 

R(t)  =  £  \  exp  [-(|)k]dx 

which  reduces  to 

R(t)  =  exp( -( t/9 )k )  (10) 

The  MLE  of  reliability  is  found  by  substituting  the  MLEs 
for  k  and  0  (k  and  0)  into  Eq  (10): 

R(t)  =  exp(-( t/0)k)  (11) 

Each  of  these  estimated  component  reliabilities  can  then  be 
used  to  estimate  a  system  reliability. 

System  Reliability 

Calculation  of  system  reliability  is  achieved  by  using 

the  laws  of  probability,  given  the  reliability  of  each 

component.  Since  it  is  assumed  that  the  components  in  each 

system  or  network  being  analyzed  fail  independently,  the 

following  equations  can  be  used  to  obtain  an  overall  system 

reliability  R  (t)  : 

s 

[Note:  Since  the  component  reliabilities  are  esti¬ 

mates,  the  system  reliability  is  also  an  estimate, 

fi8^>  1 


1.  If  the  system  is  composed  of  two  components  con¬ 
nected  in  series,  the  reliability  can  be  expressed  by 

Rs(t)  =  R1(t)-R2(t)  (12) 

A  A 

where  R^(t)  and  R2(t)  are  component  estimates. 

2.  If  the  system  is  composed  of  two  components  con¬ 
nected  in  parallel,  the  reliability  can  be  expressed  by 

Rs(t)  =  1  -  [1-R1(t)]  [1-R2(t)3  (13) 

=  1  -  QjQg  where  =  1  -  Ri(t) 

[Note:  Q(t)  represents  the  probability  that  a 

component  has  failed] 

More  complex  systems  can  be  reduced  to  combinations 
of  series  and/or  parallel  configurations  by  use  of  Bayes' 
Theorem  or  the  Boolean  Disjunctive  Theorem  (Ref  7:5-6). 

Bias  of  Reliability  Estimate 

An  important  criteria  of  an  estimator  is  its  bias.  It 
is  usually  desired  that  the  bias  be  as  small  as  possible — 
approaching  zero.  The  bias,  B  ,  is  defined  as  the  expected 
value  of  the  estimator  minus  the  true  value  of  the  parameter 
being  estimated  (Ref  11:266).  That  is 

B  =  E [R( t ) ]  -  R(t) 

In  a  Monte  Carlo  simulation  conducted  by  Thomas,  Bain,  and 
Antle  (Ref  19:365),  10,000  estimates  of  R(t)  using  sample 
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sizes  of  8  to  100  and  true  reliabilities  of  0.5  to  0.98 
were  estimated.  The  E[R(t)]  was  obtained  by  averaging 
the  10,000  estimates.  Subtracting  out  the  true  reliability 
yielded  the  bias.  Antoon  (Ref  1:44)  derived  2,000  estimates 
of  R(t)  and  calculated  the  bias  for  the  same  range  of 
sample  sizes  and  true  reliabilities  as  Thomas,  Bain,  and 
Antle.  Antoon 's  results  compared  favorably  and  are  shown 
in  Tables  I  and  II. 

Median  Rank  Values 

Given  an  ordered  random  sample  y-^ ,  y2 »  •  •  •  >  Yn  from  a 
population  having  a  Weibull  cumulative  distribution  function 
F(y )  ,  where  y  is  continuous,  estimators  can  be  deter¬ 

mined  for  F(y1),F(y2), . . . ,F(yn)  .  The  distribution  for 
these  estimators,  when  used  in  life  testing,  is  commonly 
termed  the  rank  distribution.  This  rank  distribution  is 
derived  in  the  reference  and  is  denoted  as  follows  (Ref  17: 
297-298) : 


n!  j-1  „  . 

(j-1)!  (n-j)i  pj  (l-Pj)n'JdPj  (14) 

where  0£p  .<1  and  p.  =  F(y.)  is  the  fraction  of  the 
J  J  J 

population  failing  prior  to  the  jth  ordered  observation  in 
a  sample  size  n  ,  which  by  differentiating  yields 

dF(Xj )  =  f (Xj )dx  =  dp j 
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TABLE  II 


The  median  value,  commonly  called  the  median  rank,  can 
be  found  by  considering  the  integral: 

P  -  /bCPjMPj 
0 

where  g(p.)  is  defined  in  Eq  (14).  The  value  of  x  for 
J 

which  P=0.5  would  be  the  desired  median  value.  An  approxi¬ 
mation  to  this  median  rank  value  that  was  used  in  this 
thesis  is  given  by  (Ref  17:300): 

x  =  ( j-0 . 3)/(n+0 . 4 )  (15) 


20 


III.  Procedure 


This  thesis  presents  a  modified  Double  Monte  Carlo 
method  for  obtaining  confidence  intervals  for  the  reliability 
of  various  component  networks  for  a  specific  mission  time. 
This  method  will  include  the  bias  values  and  median  rank 
derivation  given  in  Chapter  II. 

Double  Monte  Carlo  Method 

The  Double  Monte  Carlo  Method  for  obtaining  confidence 
limits  for  system  reliability  is  used,  as  a  technique,  in 
cases  where  the  distribution  (joint  distribution)  of  the 
estimator(s)  for  the  failure  model  parameter(s)  is  unknown. 

It  is  assumed  that  (Ref  13):  (1)  the  mathematical  model  for 

the  underlying  failure  distribution  is  known;  (2)  a  mathe¬ 
matical  relationship  relating  system  reliability  and  com¬ 
ponent  reliability  exists;  and  (3)  the  components  of  the 
system  have  been  subjected  to  life  tests  and  that  the  para¬ 
meters  of  the  failure  model  have  been  estimated. 

Using  the  life  test  failure  times  for  each  component, 
and  the  estimated  parameters  of  the  failure  model,  samples 
are  generated  whose  failures  have  the  same  parameters  and 
distribution,  Weibull  in  this  case,  as  the  known  or  esti¬ 
mated  ones.  It  is  important  that  these  new  samples  have  the 
same  number  of  observations  as  the  original  test  data.  The 
estimates  of  the  new  parameters  must  be  obtained  via  the 
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same  method,  maximum  likelihood,  and  with  the  same  sample 
size  used  on  the  original  sample.  These  simulated  values 
of  the  parameters  along  with  a  specified  mission  time  are 
substituted  in  the  life  distribution  to  obtain  a  simulated 
reliability,  denoted  R^  ,  for  each  component. 

In  order  to  limit  the  amount  of  computer  processing 

A 

time,  a  given  number  of  component  reliabilities,  , 

will  be  estimated  with  the  Double  Monte  Carlo  method,  for 
each  component  used  in  a  given  system.  This  given  number 
of  Ri  will  be  estimated  for  each  Monte  Carlo  iteration. 

An  empirical  distribution  will  be  formed  for  each  component 

A 

using  the  derived  component  reliabilities,  ordered 

from  smallest  to  largest,  as  the  abscissae  with  an  associ¬ 
ated  ordinate  axis  of  ordered  median  ranks.  Recall  that 
these  medium  ranks  are  calculate  using  Eq  (15)  which  is: 

x  n  +  0.4  (15) 

It  should  be  noted  that  the  first  and  last  order  statistic 
for  each  component  reliability,  associated  with  the  median 
ranks  0  and  1  respectively,  are  approximated  using  linear 
extrapolation  off  of  the  two  sequential  order  statistics 

A 

nearest  the  first  or  last  R^  .  As  an  example  of  an  empir¬ 
ical  distribution  see  Figure  2.  (Note:  Each  component 
reliability  is  determined  with  ten  sample  failures.)  These 
distributions  can  then  be  used  to  obtain  a  large  number  of 

estimated  component  reliability. 
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MEDIAN 


1.0 


Fig  2.  Sample  Empirical  Distribution  with 
Shape  =2.0  and  Scale  =  250 

Using  the  techniques  specified  in  Chapter  II,  the  true 
reliabilities  can  be  calculated  for  each  component  and  for 
each  system.  These  true  reliabilities  will  provide  an 
absolute  measure  against  which  the  modified  Double  Monte 
Carlo  method  can  be  gauged. 

Calculating/Veri fyi  ng  System  Con f idence  Limi ts 

Given  that  sample  values  of  system  reliability  have 
been  generated,  they  are  then  ordered  to  yield  a  sample 
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cumulative  distribution.  Using  this  distribution,  the 
confidence  interval  and  limits  can  be  approximated  at  any 
level  of  confidence.  For  example,  in  an  ordered  sample  of 
1,000  reliability  points,  the  one  hundredth  value  represents 
the  lower  limit  of  a  one  sided  confidence  interval  at  the 
90%  confidence  level. 

The  steps  for  finding  confidence  intervals  about  the 
true  system  reliability,  using  the  modified  Double  Monte 
Carlo  method  and  verifying  the  accuracy  of  the  confidence 
levels  are  as  follows: 

1.  Using  the  true  shape,  scale  and  location  para¬ 
meters,  generate  a  simulated  sample  of  test  data  (component 
failure  times)  from  the  Weibull  distribution. 

2.  Based  on  component  test  data,  calculate  the  maxi¬ 
mum  likelihood  estimators  of  the  shape  and  scale  parameters. 
The  location  parameter  is  assumed  to  be  zero. 

3.  Samples  are  generated  whose  failures  have  the 
same  distribution  and  parameters  as  the  estimated  parameters 
in  (2)  and  the  same  number  of  observations  as  the  original 
test  data. 

4.  The  parameters  are  again  estimated  from  the  simu¬ 
lated  sample  by  the  same  method  as  used  on  the  original 
sample  (i.e.,  maximum  likelihood  estimates). 

5.  Substitute  these  simulated  values  of  the  parameters 


into  the  reliability  function  to  obtain  a  maximum  likelihood 
estimator  of  the  component  reliability. 


6.  Subtract  the  bias  from  the  maximum  likelihood 


estimator  of  the  reliability,  using  Table  I,  to  obtain  an 
unbiased  estimator  of  the  component  reliability. 

7.  Repeat  steps  3-6  to  obtain  the  desired  number  of 
component  reliabilities  for  the  empirical  distribution. 

8.  Sample  from  this  empirical  distribution  a  given 
number  of  uniform  random  deviates  to  obtain  the  desired 
number  of  component  reliability  estimates. 

9.  Repeat  steps  1-8  for  each  component. 

10.  Calculate  point  samples  of  system  reliabilities 
from  the  point  samples  of  component  reliabilities. 

11.  Order  the  point  samples  of  system  reliabilities 
and  determine  the  99,  95,  90,  80,  70,  60,  50  percent  one¬ 
sided  confidence  intervals.  Note  if  each  of  these  intervals 
contains  the  true  system  reliability. 

12.  Repeat  step  1-11  until  the  desired  Monte  Carlo 
size  is  reached. 

13.  To  measure  the  accuracy  of  the  confidence  limits, 
determine  the  percentage  of  the  runs  in  which  each  of  the 
confidence  intervals  covered  the  true  system  reliability. 

Elaboration  on  each  individual  step  of  the  technique 
is  in  Appendix  A.  A  flow  diagram  and  computer  program  which 
executes  the  above  modified  Double  Monte  Carlo  method  are 
shown  in  Appendices  B  and  C  respectively. 


p 


Components 


Five  components  were  considered  in  testing  the  modi¬ 
fied  Double  Monte  Carlo  method.  The  true  reliability  of 


each  component,  FL  ,  is  found  by  substituting  into  the 
Weibull  reliability  formula,  Eq  (11).  Mission  time,  t  , 
is  arbitrarily  set  at  100  hours.  The  following  components 


were  used  (Ref  15:22-23): 

Component  1 

Failure  Distribution 
Parameter  Values 
True  Reliability 

Component  2 

Failure  Distribution 
Parameter  Values 
True  Reliability 

Component  3 

Failure  Distribution 
Parameter  Values 
True  Reliability 

Component  4 

Failure  Distribution 
Parameter  Values 
True  Reliability 

Component  5 

Failure  Distribution 
Parameter  Values 
True  Reliability 


Weibull 

k  =  2  0  =  250  c  =  0 

Rx  =  exp [-( 100/250) 2 ]  =  .85214 


Weibull 

k  =  3  0  =  210  c  =  0 

R2  =  exp  [-(100/210)3]  =  .89765 


Weibull 

k  =  2  6=  300  c  =  0 

R3  =  exp [-(100/300)2]  =  .89484 


Weibull 

k  =  3.5  0  =  150  c  =  0 

R4  =  exp  [-(  100/150)  3,5  ]  =  .78512 


Weibull 

k  =  2.5  0  =  250  c  =  0 

R5  =  exp[-(100/250)2'5]  =  .90376 


Systems  to  be  Analyzed 

These  five  components  were  combined  to  form  various 
kinds  of  systems.  Four  different  systems  were  used  and  are 
shown  in  Figure  3. 
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System  1 


System  3 


Pig  3.  Systems  1,  2,  3,  and  4  (Ref  15:24) 


System  1  consists  of  three  components  in  series.  The 


true  reliability  of  System  1,  R  ,  is 

S1 


Rs1  R1R2R3 


(13) 


=  (. 85214 )(. 89765) ( .89484) 
=  .68448 


System  2  consists  of  one  component  connected  in  series 
with  two  in  parallel.  Let  Qi  be  the  probability  of 
failure  for  component  i  ,  then  =  1  -  FL  .  Therefore, 

Rg  =  R1(1-Q2Q3)  (14) 

2 

=  ( .85214)  [l-( .10235)( . 10516)] 

=  .84297 


System  3  consists  of  three  components  connected  in 
parallel.  Hence, 


l  -  q1q2q3 


(15) 


=  1  -  (.14786)(.10235)(. 10516) 
=  .99841 


System  4  is  a  larger  complex  network  of  the  previous 
systems,  hence 


Rg  =  R1[l-Q2[l-R5(l-Q3Q4)n  (16) 

t>4 

=  .8521411-. 10235 [ 1 -. 90376 ( 1 - ( .10516)( .21488)))] 

=  .84197 
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IV.  Results 


It  is  known  that  the  error  associated  with  a  Monte 
Carlo  calculation  is  proportional  to  1  / »/N  where  N  is 
the  number  simulations  (Ref  17:259).  In  this  case,  the 
error  is  statistical,  that  is,  the  probable  error  is  pro¬ 
portional  to  1 / JW~  or  the  probability  is  high  that  the 
approximate  solution  does  not  deviate  from  the  true  solu¬ 
tion  by  more  than  a  certain  amount  (Ref  5:22).  The  amount 
of  statistical  error  in  a  calculation  should  decrease  with 
the  use  of  a  high  speed  digital  computer,  but  the  computer 
can  add  random  errors  due  to  arithmetic  roundoff. 

Tables  III  through  VI  contain  the  results  of  the 
modified  Double  Monte  Carlo  simulation  for  Systems  1  through 
4.  The  simulation  was  run  for  each  system  with  component 
test  data--sample  size  10,  15,  20  and  30.  In  all  cases 
run,  the  number  of  system  reliability  estimates  used  to 
partition  the  interval  [0,1]  was  1000.  The  empirical 
distribution  for  each  component  reliability  was  75  points. 
The  numbers  specified  at  each  table  entry  are  the  percent¬ 
ages  of  the  Monte  Carlo  runs  in  which  the  true  system 
reliability  was  contained  within  the  simulated  confidence 
interval.  For  example,  the  .9781  in  Table  II i  row  1  sample 
size  10,  implies  that  in  approximately  97  to  98  percent  of 
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TABLE  III 


System  1  Confidence  Interval  Coverage  of 
the  True  System  Reliability 


Confidence  Interval 
(Percent ) 

10 

Sample 

15 

Size 

20 

30 

99 

.9781 

.9831 

.9910 

.9922 

95 

.9203 

.9372 

.9371 

.9429 

90 

.8771 

.  9699 

.  8702 

.  8797 

80 

.7454 

.7667 

.7689 

.7921 

70 

.6092 

.6317 

.6209 

.6441 

60 

.4666 

.  5010 

.  5011 

.5171 

50 

.3891 

.4237 

.  3788 

.  3626 

TABLE 

IV 

System 

2  Confidence  Interval  Coverage  of 
the  True  System  Reliability 

Confidence  Interval 

(Percent)  10 

Sample 

15 

Size 

20 

30 

99 

.9600 

.9592 

.9709 

.9777 

95 

.8781 

.9000 

.9213 

.9346 

90 

.  8003 

.  8419 

.8671 

.8779 

80 

.6811 

.  7003 

.  7118 

.  7135 

70 

.5602 

.5617 

.6177 

.6389 

60 

.4765 

.  4883 

.4890 

.5049 

50 

.3617 

.  3812 

.  3900 

.4371 

TABLE  V 


System  3  Confidence  Interval  Coverage  of 
the  True  System  Reliability 


Confidence  Interval 
(Percent ) 

10 

Sample 

15 

Size 

20 

30 

99 

.8598 

.9017 

.9221 

.9486 

95 

.7001 

.  7786 

.7983 

.8286 

90 

.5410 

.6443 

.664  7 

.7057 

80 

.  3920 

.4892 

.4701 

.5577 

70 

.2871 

.  3337 

.  3686 

.4171 

60 

.1999 

.2209 

.2300 

.3057 

50 

.  1321 

.14  37 

.1611 

.2251 

TABLE  VI 

System  4  Confidence  Interval  Coverage  of 
the  True  System  Reliability 


Confidence  Interval 
(Percent ) 

10 

Sample 

15 

Size 

20 

30 

99 

.9553 

.9590 

.9692 

.9735 

35 

.8807 

.9011 

.9273 

.9440 

90 

.8176 

.  8348 

.8600 

.8900 

80 

.7186 

.7019 

.7231 

.7373 

70 

.5719 

.5882 

.6147 

.6306 

60 

.4662 

.4760 

.4081 

.4779 

50 

.3797 

.4010 

.  3980 

.4207 

the  runs  for  System  1,  the  99  percent  confidence  interval 
contained  the  true  reliability.  It  can  be  observed  that, 
within  the  limits  of  the  simulation  error,  as  the  sample 
size  increases  for  all  the  systems,  the  confidence  interval 
coverage  improves.  In  those  few  cases  where  it  did  not 
improve,  variability  of  the  method  can  be  seen. 

There  is  a  noticeable  effect  of  the  true  system  relia¬ 
bility  on  the  confidence  interval  coverage.  It  appears 
that  a  lower  system  reliability  corresponds  to  a  more  accu¬ 
rate  confidence  interval.  System  1,  with  a  reliability  of 
.68448,  has  consistently  more  accurate  confidence  interval 
coverage  than  any  other  system.  System  2  and  4,  with 
reliabilities  of  .84297  and  .84197  respectively,  have  table 
values  that  correspond  closely.  System  3,  with  the  highest 
reliability,  .99841,  has  the  least  accuracy. 

As  an  indication  of  the  time  that  was  required  to  run 
these  simulations  on  the  CDC  Cyber,  Table  VII  specifies 
for  all  systems  together,  at  different  sample  sizes,  the 
average  amount  of  time  12,  50  iteration  program  runs,  to 
obtain  600  Monte  Carlo  iterations.  That  is,  each  program 
that  was  run  contained  all  4  systems  at  a  given  sample  size, 
for  50  iterations. 

Tables  VIII  through  XI  provide  a  comparison  of  the 
Modified  Double  Monte  Carlo  method  with  the:  Univariate 
Method  (Ref  15:33-34)  for  all  systems  at  sample  size  20; 
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TABLE  VII 


CPU  Times  on  the  CDC  Cyber 


Sample 

Size 

10 

15 

20 

30 

1456  sec 

1871  sec 

2378  sec 

2982  sec 

Bivariate  Method  (Ref  5:29)  System  4  (only)  sample  size  20 
for  700  Monte  Carlo  iterations. 

For  sample  size  20,  the  modified  Double  Monte  Carlo 
provides  a  more  consistent  coverage  of  the  system  relia¬ 
bility  at  confidence  levels  of  99,  95  and  90.  At  the 
remaining  confidence  levels,  the  univariate  is  more  accurate. 
A  reason  for  this  is  that  since  the  distribution  of  the 
estimators  for  the  parameters  is  unknown  in  the  Double 
Monte  Carlo,  at  the  lower  confidence  limits  the  distributioh 
is  such  that  more  reliability  points  appears  in  the  distri¬ 
bution  tails.  It  can  also  be  noted  that  in  most  cases 
the  confidence  interval  coverage  is  low  for  small  sample 
sizes.  Since  there  is  high  or  optimistic  system  reliability 
estimates,  the  lower  confidence  limits  are  also  too  high 
and  the  confidence  interval  is  not  wide  enough  to  cover  the 
true  system  reliability. 

Comparing  the  modified  Double  Monte  Carlo  method  with 
the  bivariate  method,  it  can  be  seen  that  the  bivariate 


TABLE  VIII 


System  1  Comparison  of  Modified  Double 
Monte  Carlo  and  Univariate 


Confidence  Interval 
(Percent ) 

Univariate 

Double  Monte  Carlo 

99 

.9617 

.9910 

95 

.9050 

.9371 

90 

.  8617 

.8702 

80 

.  7800 

.  7689 

70 

.6817 

.6709 

60 

.5700 

.5011 

50 

.4833 

.3788 

TABLE  IX 

System  2  Comparison  of  Modified  Double 
Monte  Carlo  and  Univariate 


Confidence  Interval 
(Percent ) 

Univariate 

Double  Monte  Carlo 

99 

.9383 

.9709 

95 

.8600 

.9213 

90 

.8083 

.8671 

80 

.7200 

.7118 

70 

.6467 

.6177 

60 

.  5433 

.4890 

50 

.4700 

.  3900 
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TABLE  X 


System  3  Comparison  of  Modified  Double 
Monte  Carlo  and  Univariate 


Confidence  Interval 
(Percent ) 

Univariate 

Double  Monte  Carlo 

99 

.9017 

.9221 

95 

.7950 

.7983 

90 

.7050 

.6647 

80 

.5733 

.4701 

70 

.4617 

.3686 

60 

.  3650 

.2300 

50 

.2767 

.1611 

TABLE  XI 

System  4  Comparison  of  Modified 
with  Univariate  and 

.  Double  Monte 
Bivariate 

Carlo 

Confidence  Interval 
(Percent ) 

Univariate 

Double 

Monte  Carlo 

Bivariate 

99 

.9417 

.9692 

.996 

90 

.  8033 

.  8600 

.926 

80 

.7183 

.7231 

.826 

70 

.6450 

.6147 

.730 

60 

.5533 

.4681 

.604 

50 

.4733 

.  3980 

.490 
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ethod  is  conservative,  confidence  interval  is  wider, 
thereby  providing  a  more  accurate  confidence  level  that  is 
much  less  sensitive  to  degradation  due  to  high  system 
reliability . 
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V.  Conclusions  and  Recommendations 


Conclusions 

The  modified  Double  Monte  Carlo  method  developed  in 
this  thesis  provides  a  useful  tool  for  reliability  estima¬ 
tion  when  there  is  a  lack  of  component  test  data,  that  is, 
when  there  are  limited  test  times  to  failure  for  a  given 
component.  The  modified  Double  Monte  Carlo  demonstrates 
better  results  with  respect  to  confidence  levels  for  system 
reliability  than  the  univariate  case,  when  the  sample  size 
is  small  and  a  higher  degree  of  confidence  is  required.  As 
with  the  univariate  case,  if  there  is  a  high  component  or 
system  reliability,  the  confidence  intervals  capture  the 
true  reliability  to  a  lesser  extent. 

It  can  also  be  concluded  that  the  skewness  of  the 
empirical  distribution  of  the  component  reliabilities  R(t) 
has  a  definite  impact  on  the  results,  to  the  extent  that  at 
the  lower  confidence  intervals  more  reliabilities  are  out¬ 
liers  in  the  distribution  tails. 

The  modified  Double  Monte  Carlo  method  presented  in 
this  study  will  enable  anyone  to  obtain  approximate  confi¬ 
dence  intervals  when  the  failure  models  are  from  the  one, 
two  or  three  parameter  Weibull  and  the  distribution  of  the 
estimators  for  the  parameters  are  unknown. 
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Recommendations 


The  modified  Double  Monte  Carlo  should  be  used  to 
investigate  confidence  levels  when  the  component  reliabil¬ 
ity,  R(t)  ,  failure  times  are  modeled  by  the  exponential, 
normal  or  lognormal  distribution.  Again,  this  should  be 
attempted  when  the  sample  size  is  relatively  small.  Emphasis 
could  also  be  placed  on  component  reliabilities  greater  than 
0.9. 

Because  of  the  limitation  on  the  component  empirical 
distribution,  sensitivity  analysis  should  be  conducted  to 
see  what  effect  might  take  place  when  the  distribution  size 
is  increased  to  150,  200,  and  300  points. 
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Appendix  A 
Notes 


The  following  notes  pertain  to  the  steps  outlined  in 
the  modified  Double  Monte  Carlo  procedure  in  Chapter  III. 

Step  1 :  Weibull-distribution  sample  failure  times 
were  generated  by  using  the  International  Mathematical  and 
Statistical  Libraries  (IMSL)  subroutine  GGWIB. 

Step  2:  The  Maximum  Likelihood  Estimators  (MLE)  of 

A  A 

the  shape  and  scale  parameters  (k  and  8)  were  determined 
by  use  of  a  general  solution  subroutine  of  Weibull  para¬ 
meters  called  PARES  (Ref  7).  In  this  case,  NSAM  is  the 
failure  sample  size,  m  =  NSAM  (i.e.,  there  are  no  obser¬ 
vations  remaining  after  censoring),  MR  =  0  no  observations 
are  censored,  R  is  the  I-th  order  statistic  of  sample 
(1=1,  NSAM)  .  (Note:  This  iterative  procedure,  in  the 
case  of  a  Weibull  population,  is  applicable  to  the  most 
general  case  in  which  all  three  parameters  are  unknown  and 
must  be  solved  simultaneously.  It  is  also  applicable  to 
special  cases  in  which  any  one  or  any  two  of  the  parameters 
are  unknown.  This  is  accomplished  by  specifying  combina¬ 
tions  of  ssl  ,  ss2  ,  or  ss3  equal  to  one  or  zero.) 

Steps  3 ,  4 :  Double  Monte  Carlo  step  where  the  compo- 

-A  A 
A  A 

nent  parameters  (k,G)  are  reestimated  via  MLE. 

Step  5:  The  MLE  of  component  reliability,  R(t)  , 
was  found  by  substitution  of  k  and  0  into  Eq  (11). 


Step  6:  The  bias  of  the  MLE  R(t)  was  determined 
by  interpolation  (using  cubic  splines)  in  Table  I. 

Step  7:  The  first  and  last  order  statistics  of  each 
component  reliability  were  determined  by  linear  extrapola¬ 
tion  on  the  two  nearest  reliability  points.  Subroutine 
EXTRA  is  used  for  this  purpose. 

Step  8:  A  vector  of  1000  component  reliabilities  is 
obtained  from  the  empirical  distribution  using  an  IMSL 
subroutine  GGUBS  to  generate  uniform  random  deviates. 

Cubic  splines  are  used  to  interpolate  for  the  actual  com¬ 
ponent  reliabilities. 

Step  10:  The  vectors  of  component  reliabilities  are 
combined  using  the  system  reliability  equations  [Eqs  (13), 
(14),  (15),  and  (16)].  This  results  in  a  vector  of  1000 
estimated  system  reliabilities  for  each  system. 

Step  11:  These  system  reliability  vectors  are 
ordered  in  ascending  sequence  to  obtain  the  1,  5,  10,  20, 

30,  40,  50  percentiles. 

Step  13:  If  the  lower  confidence  limit  is  less  than 
or  equal  to  the  true  system  reliability  then  the  interval 
contains  the  true  reliability,  subroutine  C0NLIM  is  used 
for  this  purpose. 

In  this  thesis,  500  Monte  Carlo  runs  of  the  simulation 
were  made  in  Step  12.  For  each  run  and  system,  it  was  noted 
if  the  true  reliability  was  contained  in  the  confidence 
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intervals.  As  a  comparison,  the  X  percent  confidence 
interval  should  contain  the  true  system  reliability  X 
percent  of  the  time.  This  is  the  validation  of  the  method. 
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Appendix  C 

Computer  Program  Listing 


PROGRAM  DOUMC  (INPUT,  OUTPUT) 

C 

c  ****************************************************************** 
C  *  THIS  PROGRAM  GENERATES,  FOR  EACH  COMPONENT  OF  ANY  SYSTEM,  A  * 
C  *  SAMPLE  SIZE  NS AM  OF  FAILURE  TIMES  FROM  THE  WE I  BULL  DISTRIBUTION* 
C  *  USING  THE  TRUE  PARAMETERS  TK( SHAPE),  TTHETA( SCALE )  AND  TC  * 

C  *  (LOCATION).  FROM  THIS  SAMPLE,  THE  MAXIMUM  LIKELIHOOD  * 

C  *  ESTIMATORS  ( MLES )  FOR  K  AND  THETA  ARE  DERIVED  USING  HARTER  AND  * 
C  *  MOORE  ITERATIVE  SCHEME ( PARES ) .  THESE  ESTIMATES  FOR  K  AND  THETA* 
C  *  ARE  THEN  USED  TO  GENERATE  A  NEW  SAMPLE  SIZE  NS AM  OF  WIKBULL  * 
C  *  FAILURES,  WHICH  ARE  IN  TURN  USED  TO  RECALCULATE  THE  MLES  FOR  K  * 
C  *  AND  THETA.  THESE  SECOND  MLES  ARE  COMBINED  TO  YIELD  CHAT,  THE  * 
C  *  MLE  FOR  THE  COMPONENT  RELIABILITY.  BIASES  FOR  THESE  ESTIMATES  * 
C  *  ARE  OBTAINED  AND  SUBTRACTED.  A  GIVEN  NUMBER  (DSAM)  OF  THESE  * 
C  *  UNBIASED  CHAT  ARE  ESTIMATED  WHICH  ARE  USED  TO  FORM  THE  X-AXIS  * 
C  *  FOR  AN  IMIPKICAL  DISTRIBUTION.  THE  Y-AXIS  IS  OBTAINED  USING  A  * 
C  *  MEDIAN  RANK  PROCEDURE.  THIS  E.MIPR1CAL  DISTRIBUTION  IS  THAN  * 
C  *  USED  TO  SAMPLE  (NGE.N)  NUMBER  OF  COMPONENT  RELIABILITIES.  THIS  * 
C  *  PROCESS  IS  REPEATED  FOR  EACH  COMPONENT  AND  THEN  RELIABILITIES  * 
C  *  ARE  COMBINED  FOR  4  DIFFERENT  SYSTEMS  TO  YIELD  4  VECTORS  OF  * 

C  *  SAMPLE  SYSTEM  RELIABILITIES.  THESE  VECTORS  ARE  THEN  ORDERED  * 
C  *  AND  THE  99, 95,  90, 80, 70, 60,  SO  PERCENT  LOWER  CONFIDENCE  LIMITS  * 
C  *  ARE  PICKED  (CONLIM).  IT  IS  NOTED  IF  THE  TRUE  RELIABILITY  IS  * 
C  *  CONTAINED  WITHIN  THESE  INTERVALS.  * 


C  *  THE  ABOVE  PROCESS  IS  REPEATED  FOR  NOI.MC  MONTE  CARLO  RUNS  WITH  * 
C  *  COUNTERS  FOR  EACH  SYSTEM  TO  TRACK  THE  NUMBER  OF  TIMES  THAT  THE  * 
C  *  TRUE  RELIABILITY  IS  CAPTURED.  * 

c  ****************************************************************** 
DIMENTION  ARC  1 1  AT ( 1 000 , 5 ) , ATRS ( 4 ) , B I ( 12,15) ,BPAR(4) ,C( 4,4 ) , 
1CCHAT ( 100) ,PARAM(3,5 ) ,R( 500,5) ,RANK( 100) ,RI.BS( IB ) ,SCHAT( 1000) , 
2SMSZ( 15) , SPLINE ( 99,3) ,TEMP( 100) ,TRC(5) ,TRS( 1 ) ,UNIF( 1000) , 

3WK( 200 ) ,X( 2 ) , ZCHAT ( 100) ,CC( 550) , THETA ( 550) ,EK( 550) 

DOUBLE  PREC ISIONDSEED 
INTEGER  DSAM 
DATA  BPAR/O . , 0 . , 0 . , 0 . / 

DATA  RLBS/.5 ,  .  55  ,  . 6  ,  . 55 , . 7 , . 75 , . 8 , . 85 , . 9 , . 925 , . 95 , 

1  .98/ 

DATA  SMSZ/8 . , 9 . , 1 0 . , 1 1 . , 1 2 . , 1 3 . , 1 4 . , 1 5 . , 20 . , 25 . , 

1  30. ,40. ,50. ,75. , 100./ 

C 

C  INITIALIZE  VARIABLES*** 

READ* .DSEED , T 
CALL  RANSET(T) 

TIME  -100. 

NRLBS-12 
NSMSZ-- 1  5 
NMC 1-2 

NSMSZ 1 =NSMSZ- 1 
NRI.BS1  -  NRI.BS- 1 
ISD  12 
NGEN- 1000 
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NS1C99  =  0 
NS1C95  =  0 
NS1C90  =  0 
NS1C80  =  0 
NS1C70  =  0 
NS1C60  =  0 
NS1C50  =  0 
NS2C99  =  0 
NS2C95  =  0 
NS2C90  =  0 
NS2C80  =  0 
NS2C70  =  0 
NS2C60  =  0 
NS2C50  =  0 
NS3C99  =  0 
NS3C95  =  0 
NS3C90  =  0 
NS3C80  =  0 
NS3C70  =  0 
NS3C60  =  0 
NS3C50  =  0 
NS4C99  =  0 
NS4C95  =  0 
NS4C90  =  0 
NS4C80  =  0 
NS4C70  =  0 
NS4C60  =  0 
NS4C50  =  0 
C 

WRITE  400 
C 

C  READ  SAMPLE  SIZE  AND  NUMBER  0E  MONTE  CARLO  RUNS*** 

C 

READ  * , NSAM , NOLMC , DSAM 

WRITE  401 , NSAM, NOLMC , DSAM 
C 

C  SET  PARAMETERS  FOR  PARES  AND  RANK  DISTRIBUTION*** 

C 

M=NSAM 

MR=0 

RNSAM- NSAM 

NNMC2 -DSAM-2 

TSAM=  DSAM 

IC-DSAM— 1 
C 

C  READ  TRUE  COMPONENT  PARAMETERS*** 

C 

READ  *,  ( ( PARAM  (  I  ,  .1 )  ,  I  -- 1 , 3  )  ,  J  - 1  ,  S  ) 

WRITE  402,  (  (J,  (  PARAM  (  I  ,  J  )  ,  I  1  ,3  )  ,CMPREI.(  TIME  ,  PARAM  (  1  ,.I)  , 
1PARAMI 2 , J ) , PARAM ( 3 , J ) ) ) ,J+1 ,5) 

C 

C  READ  BIAS  LEVELS*** 

C 
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READ  * ,  ( ( B I ( I ,J ) , J - 1 ,15) ,1  =  1 , 12) 
WRITE  500 

WRITE  403,((BI(I,J),J=1,15) ,1=1,12) 

CALCULATE  MEDIAN  RANK  DISTRIBUTION*** 
C 

RANK( 1 ) -0 . 

RANK(DSAM)  1 . 

DO  305  1=1 , NNMC2 

RANK( I + 1 ) = ( I - . 3 ) / ( TSAM+ . 4 ) 

305  CONTINUE 


C  THIS  LOOP  (FROM  HERE  TO  STATEMENT  300)  COMPLETES  NOLMC  MONTE  CARLO 
C  RUNS  OK  THE  PROGRAM*** 

C 


DO  300  NCOUNT  1, NOLMC 
C 

C**************************'****  *•**★**■*■*■****★*  vV  A********  *************** 

C  FOR  EACH  OF  5  COMPONENTS,  THIS  LOOP  GENERATES  THE  MI.E 
C  OF  CHAT  FROM  THE  MI.ES  OF  K  AND  THETA.  IT  SUBTRACTS  THE 
C  BIAS:  ESTABLISHES  THE  EM  1  PR  I  CAL  DISTRIBUTION  AND  SAMPLES 
C  THE  GIVEN  NUMBER  OF  COMPONENT  RELIABILITIES*** 

DO  200  J  1,5 
TK  PARAM (  1  ,.)) 

TC  PARAM ( 3 , J ) 

C 

C  DETERMINE  TRUE  COMPONENT  RELIABILITY*** 

TRC(J)  CMPREI.  ( T I  ME ,  TK ,  TTHET*  C) 

C 

C  THIS  LOOP  GENERATES  THE  DESIRED  NUMBER  OF  EMIPRICAL 
C  DISTRIBUTION  POINTS*** 

C 

DO  150  LI.  1.NNMC2 
LK=  LI. 

C 

C  GATHER  NS AM  FAILURE  TIMES  FROM  THE  KE I  BULL 
C  DISTRIBUTION  WITH  TRUE  PARAMETERS  TK , TTHETA 
C  AND  TC*** 

220  CALL  GGWI B ( DSEED , TK , NS AM , TEMP ) 

DO  201  I  1 , NSAM 
R( I , J)  TTHETA*TEMP( I ) tTC 
201  CONTINUE 
CC< 1 ) :  TC 
THETA ( 1 )  TTHETA 
EK( 1 )  TK 
C 

C  DETERMINE  THE  MI.ES  FOR  THETA  AND  K*** 

CALL  PARES(NSAM,M,CC, THETA, EK, MR, TTHETA, TK,R( I , J ) ) 

IFILK.NE. 1 )G0  TO  215 
X( 1 ) = TTHETA 
X(2)-TK 
LK=LK+ 1 
GO  TO  220 
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c 

C  DETERMINE  THE  COMPONENT  RELIABILITY  USING  THE  TRUE  C 
C  AND  THE  TWICE  MLES  OF  K  AND  THETA*** 

215  CHAT=CMPREL (TIME , TK , TTHETA , TC ) 

TTHETA=X( 1) 

TK=X( 2 ) 

C 

C  GIVEN  THE  MLE  OF  THE  COMPONENT  RELIABILITY  AND  THE 
C  SAMPLE  SIZE,  ENTER  THE  2-DIMENSIONAL  ARRAY  BI 
C  AND  FIND  THE  BIAS  OF  THE  ESTIMATOR*** 

DO  203  1=1 , NSMSZ 1 
LSMSZ=1 

IF ( RNSAM . LE . SMSZ ( I +1 ) ) GO  TO  204 

203  CONTINUE 
GO  TO  205 

204  1 F ( CHAT . LT . . 5 ) GO  TO  206 
DO  207  1=1 ,NRLBS1 
LRLBS+I 

IF ( CHAT . LE . RLBS ( I +1 ) )G0  TO  208 

207  CONTINUE 
GO  TO  205 

208  CALL  IBCICU ( Bl , I SD , RLBS , NRLBS , SMSZ , NSMSZ , LRLBS , LSMSZ , C , WK , I ER ) 
CALL  IBCEVU( RLBS, NRLBS, SMSZ, NSMSZ, LRLBS, LSMSZ, C, CHAT, RNSAM, 

1BIAS.IER) 

I F ( LRLBS . LT .NRLBS ) GO  TO  209 

205  BIAS=0 . 

GO  TO  209 

206  IF (RNSAM. EQ. SMSZ ( LSMSZ +1 ) H.SMSZ=LSMSZ+1 
BI AS=BI ( 1 , LSMSZ ) 

C 

C  SUBTRACT  THE  BIAS*** 

209  ZCHAT(LL)=CHAT-BIAS 
150  CONTINUE 

C 

C  ORDER  THE  RELIABILITIES  FOR  THE  EMIPRICAL  DISTRIBUTION*** 

CALL  VSRTA ( ZCHAT , NNMC 2 ) 

DO  211  II, NNMC 2 
211  CCHATl I + 1 )  ZCHAT ( I  ) 

C 

C  FIND  THE  FIRST  AND  LA  IT  ORDER  STATISTICS  FOR  THE 
C  EMIPRICAL  DISTRIBUTION  COMPONENT  RELIABILITIES*** 

CALL  EXTRA ( CCHAT , RANK , NMC 1 ,0.,T2,DSAM) 

CCH AT ( 1 )  T2 

I F ( T2 . LE . 0 . ICCHATI 1 i  0. 

CALL  EXTRA <  CCILYi  , RANK, NNMC 2, 1 . ,T2 , DSAM ) 

CCHAT(DSVM)  T2 

IF(T2.GE. 1 . ICCHATI DSAM I  I . 

C 

C  SAMPLE  FROM  A  UNIFORM  GENERATOR  INGEN)  NUMBER  OF 
C  COMPONENT  RELIABILITIES*** 

CALL  GGUBS I DSEED , NGEN , I N I F  ) 

CALL  ICSICIK  RANK, CCHAT, DSAM, BPAR, SPLINE, I C, I ER ) 

CALL  ICSEVUI RANK, CCHAT , DSAM , SPL I NE , I C , UN  I F , SCIIAT , NGEN , I ER I 
DO  210  KK  I, NGEN 
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210  ARCHAT ( KK, J ) =SCHAT ( KK) 

200  CONTINUE 
C 

C  FOR  EACH  OF  THE  4  SYSTEMS,  THE  DIFFERENT  COMPONENTS  ARE 
C  COMBINED  TO  YIELD(NGEN)  SAMPLES  OF  THE  SYSTEM  RELIABILITIES. 

C  THESE  SAMPLES  ARE  SEQUENCED  IN  ASCENDING  ORDER  THE 
C  COUNTERS  KEEP  TRACK  OF  WHEN  THE  99,95,90,80,70,60,50 
C  PERCENT  CONFIDENCE  INTERVALS  CONTAIN  THE  TRUE  SYSTEM 
C  RELIABILITIES*** 

C 

C  SYSTEM  1*** 

NMC  =  1 

CALL  REL1 (NMC ,TRC ,TRS) 

ATRS ( 1 ) =TRS ( 1 ) 

CALL  REL1 (NGEN, ARCHAT, SCHAT) 

CALL  VSRTA( SCHAT, NGEN) 

CALL  CONLIM (SCHAT,TRS,NS1C99,NS1C95,NS1C90,NS1C80,NS1C70,NS1C60, 
1NS1C50 ,NGEN ) 

C 

C  SYSTEM  2*** 

NMC  =  1 

CALL  REI.2  ( NMC  ,  TRC  ,  TRS ) 

ATRS(2)-TRS( 1 ) 

CALL  REL2 ( NGEN , ARCHAT , SCHAT ) 

CALL  VSUTAt SCHAT, NGEN) 

CALL  CONLIM (SCHAT, TRS, NS2C99,NS2C95,NS2C90,NS2C80,NS2C70,NS2C60, 
1NS2C50.NGEN) 

C 

C  SYSTEM  3*** 

NMC  =  1 

CALL  REI.3  ( NMC  ,  TRC  ,  TRS ) 

ATRS( 3 ) -TRS ( 1 ) 

CALL  REL3 ( NGEN , ARCHAT , SCHAT ) 

CALL  VSRTA( SCHAT, NGEN) 

CALL  CONLIM (SCHAT, TRS, NS3C99 , NS3C95 ,NS3C90 ,NS3C80 , NS3C70 , NS3C60 , 
1NS3C50 ,NGEN ) 

C 

C  SYSTEM  4*** 

NMC=1 

CALL  REL4 ( NMC , TRC , TRS ) 

ATRC ( 4 ) =TRS ( 1 ) 

CALL  REL4 ( NGEN , ARCHAT , SCHAT ) 

CALL  VSRTA( SCHAT, NGEN) 

CALL  CONI, IM ( SCHAT , TRS , NS4C 99 , NS4C 95 , NS4C90 , NS4C80 , NS4C70 , NS4C60 , 
1NS4C50.NGEN) 

300  CONTINUE 
C 

RNOI.MC  =  NOLMC 

C  FOR  SYSTEM  1,  DETERMINE  THE  95,  90,  AND  80  PERCENT  CONFIDENCE 
C  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 
PS  1C 99  --  NS1C99/RN0LMC 
PS1C95  =  '  >  1C 95 /RNOI.MC 
PS  1C  90  =  NS1  C90/RN0I.MC 
PS  1C 80  -  NS1C80/RN0LMC 
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PS1C70  =  NS1C70/RN0I.MC 
PS1C60  =  NS1C60/RN0LMC 
PS1C50  =  NS1C50/RN0LMC 

C  FOR  SYSTEM  2,  DETERMINE  THE  95,  90,  AND  80  PERCENT  CONFIDENCE 
C  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 
PS2C99  =  NS2C99/RN0LMC 
PS2C95  =  NS2C95/RN0LMC 
PS2C90  =  NS2C90/RN0LMC 
PS2C80  =  NS2C80/RN0LMC 
PS2C70  =  NS2C70/RN0LMC 
PS2C60  =  NS2C60/RN0LMC 
PS2C50  =  NS2C50/RN0LMC 

C  FOR  SYSTEM  3,  DETERMINE  THE  95,  90,  AND  80  PERCENT  CONFIDENCE 
C  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 
PS3C99  =  NS3C  99/RNOLMC 
PS3C95  =  NS3C95/RN0LMC 
PS3C90  =  NS3C90/RN0LMC 
PS3C80  =  NS3C80/RN0LMC 
PS3C70  =  NS3C70/RN0LMC 
PS3C60  =  NS3C60/RN0LMC 
PS3C50  =  NS3C50/RN0LMC 

C  FOR  SYSTEM  4,  DETERMINE  THE  95,  90,  AND  80  PERCENT  CONFIDENCE 
C  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 
PS4C99  =  NS4C99/RN0LMC 
PS4C95  =  NS4C95/RN0I.MC 
PS4C90  =  NS4C  90/RN0LMC 
PS4C80  =  NS4C80/RN0LMC 
PS4C70  =  NS4C70/RN01.MC 
PS4C60  =  NS4C60/RN0LMC 
PS4C50  =  NS4C50/RN0LMC 

PRINT  404 

PRINT  405 , ATRS ( 1 ) , PS1C99 ,PS 1 C95 ,PS 1C90 , PS1C80 , PS  1  CTO , PS  1 C60 , PS  1 C50 

PRINT  406 

PRINT  405 , ATRS ( 2 ) ,PS2C99 , PS2C95 , PS2C90 , PS2C80 , PS2C70 , PS2C60 , PS2C50 

PRINT  407 

PRINT  405, ATRS (3 ) , PS3C 99 , PS3C95 , PS3C90 , PS3C80 , PS2C70 , PS2C60 , PS2C50 

PRINT  408 

PRINT  405 ,ATRS(4) ,PS4C99 , PS4C95 , PS4C90 , PS4C80 , PS4C70 , PS4C60 , PS4C50 

STOP 

C 

£*********★**********★********** ************************************** 

£  ★a******************  *  ************  *  *****  ‘.V  *  *  *  *  *  ******  *  *  *  *  *  *  ********  Vf  *  *  * 

C400  FORMAT* "1") 

401  FORMAT  (  "  ***************************************************" 


1 

"**********" 

,  / 

t » 

* 

*",  T62 ,  ’ 

/  , 

M  **» 

1 

T62 ,  "*", 

/ 

2 

,  "  T25 , 

"SAMPLE 

"  "SIZE 

",  13, 

T62  , 

"*",  /  .  " 

*»» 

3 

T62 ,  / 

1  1 

1 

* 

T22 ,  "MONTE  CARLO 

"  "SIZE  ",  13, 

4 

T62, 

I 

,  TIM 

,  "DISTRIBUTION  S 

IZE 

"  7  7 

t  *  '*  » 

5 

T62 ,  / 

9 

T62 , 

/  ,  "  T62 

»»*"  I 

>  »  /  9 

7  ///////) 

402  FORMAT  (5( IX,  "COMPONENT  ",  II,  /  ,  6X,  "K  ",  F4.2,  /  , 

1  6X ,  "THETA  ",  F 5.0,  /  ,  6X,  "C  ",  F2.0,  /  ,  6X , 

2  "RELIABILITY  ",  F7.5,  /  /  )  ) 
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403  FORMAT  (  15(1X,  F7.5)  ) 

404  FORMAT  (  /  "  *****  SYSTEM  1  *****"  /  "  (3  COMPONENTS  IN  SERI" 

1  "ES)"  ) 

405  FORMAT  (  /  ,  "  TRUE  SYSTEM  RELIABILITY  =" ,  F7.5,  / 

1  "  THE  99  PERCENT  CONFIDENCE  INTERVAL  COVERED  ",  F6 . 4 , 

2  "  Of  THE  RUNS",  /  ,  "  THE  95  PERCENT  CONFIDENCE  INTERVAL  COV" 

3  "ERED  ",  F6.4  "OFTHERUNS",  /  ,  "  THE  90  PERCENT  CONFIDENCE  I" 

4  "NTERVAL  COVERED  "  F6.4,  "  OF  THE  RUNS",  /  ,  "  THE  80  PERCE" 

5  "NT  CONFIDENCE  INTERVAL"  "COVERED  "  F6.4  "OFTHERUNDS" ,  /  , 

6  "  THE  70  PERCENT  CONFIDENCE  "  "INTERVAL  COVERED  ",  F6.4, 

7  "  OF  THE  RUNS",  /  ,  "  THE  60  PERCENT  "  "CONFIDENCE  INTERVAL" 

8  "  COVERED  ",  F 6 . 4 ,  "  OF  THE  RUNDS" ,  /  ,  "  THE",  "  50  PERCEN" 

9  "T  CONFIDENCE  INTERVAL  COVERED  ",  F6.4,  "  OF  THE  RUNS",  /  / 

9  /  ) 

406  FORMAT  (  /  "  *****  SYSTEM  2  *****"  /  "  (1  COMPONENT  IN  SERIE" 

1  "S  WITH  2  "  "IN  PARALLEL)"  ) 

407  FORMAT  (  /  "  *****  SYSTEM  3  *****"  /  "  (3  COMPONENTS  IN  PARA" 

1  "LLEL)"  ) 

408  FORMAT  (  /  "  *****  SYSTEM  4  *****"  /  "  (A  5-COMPONENT  COMPLE" 

1  "X  NETWORK)"  ) 

500  FORMAT  ("  ",  "BIAS  OF  ESTIMATED  COMPONENT  RELIABILITY") 

END 
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FUNCTION  CMPREL ( T IME , K , THETA , C ) 

C 

C  THIS  FUNCTION  ROUTINE  CALCULATES  THE  COMPONENT 
C  RELIABILITIES  FROM  THE  WE I BULL  DISTRIBUTION*** 
REAL  K 

CMPREL=EXP ( - ( ( T IME-C ) /THETA ) **K ) 

RETURN 

END 


SUBROUTINE  REL1 ( NMC , R , RS ) 

C  REL1  DETERMINES  THE  SYSTEM  RELIABILITY  OF  3  COMPONENTS  IN  SERIES*** 
DIMENSION  R ( NMC , 5 ) ,  RS9NMC ) 

DO  10  1=1, NMC 

10  RS( I ) =R( I , 1 )  *  R ( I , 2 )  *  R{ I , 3 ) 

RETURN 

END 


SUBROUTINE  REI.2  (  NMC  ,  R  ,  RS ) 

C  REL2  DETERMINES  THE  SYSTEM  RELIABILITY  OF  I  COMPONENT  IN  SERIES 
C  WITH  2  IN  PARALLEL*** 

DIMENSION  R ( NMC , 5 ) ,  RS(NMC) 

DO  10  1=1 ,NMC 

10  RS(I)=R(I ,1)*( 1 .-( 1 . — R ( I , 2 ) )*( 1 ,-R(I ,3) ) ) 

RETURN 

END 


SUBROUTINE  REI.3  ( NMC  , R  ,  RS  ) 

C  RE  1, 3  DETERMINES  THE  SYSTEM  RELIABILITY  OF  3  COMPONENTS  IN  PARALLEL*** 
DIMENSION  R ( NMC , 5 )  ,  RS(NMC) 

DO  10  1=1, NMC 

10  RS( I )  =  1 .-( 1 . -R ( I , 1 ) )*( 1 .-R( 1 ,2) )*( 1 .-R( I  ,3)  ) 

RETURN 

END 


54 


SUBROUTINE  REL4 ( NMC , R , RS ) 

C  REL4  DETERMINES  THE  SYSTEM  RELIABILITY  OF  A  5  COMPONENT  COMPLEX 
C  NETWORK*** 

DIMENSION  R ( NMC , 5 ) ,  RS(NMC) 

DO  10  1=1, NMC 

RS(I)=R(I,1)*(1.-(1.-R(I,2) )*(1.-R(I,5)*(1.-(1.-R(I,3) ) 
1*(1.-R(I,4) ) ) ) ) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  EXTRA( X , Y , I ,T1 ,T2 , NNMC ) 

C 

C  THIS  SUBROUTINE  IS  USED  TO  APPROXIMATE  THE  FIRST  AND  LAST  ORDER 
C  STATISTICS  OF  THE  COMPONENT  RELIABILITY  FOR  THE  EMIPRICAL 
C  DISTRIBUTION*** 

DIMENSION  X ( NNMC ) , Y ( NNMC ) 

SLOPE= ( Y ( I + 1 ) -Y ( I ) ) / ( X ( I + 1 ) -X ( I ) ) 

B=Y ( I ) -SLOPE*X( I ) 

T2= ( T1 -B ) /SLOPE 

RETURN 

END 


C 

C 

C 


SUBROUTINE  C0NLIM(R,RS,N1 ,N2 ,N3 ,N4,N5 ,N6 ,N,7 ,NMC ) 


THIS  SUBROUTINE  KEEPS  TRACK  OF  THE  NUMBER 
RELIABILITIES  IS  CONTAINED  WITHIN  A  GIVEN 
DIMENSION  R(NMC) ,RS( 1 ) 

IF( R(  10 )  .LF,.  RS(ll)  NUN1  +  1 


IF ( R( 50) 

I F ( R ( 100) 
IF(R(200) 
IF ( R ( 300 ) 
IF( R ( 400 ) 
IF(R( 500) 
RETURN 
END 


.LE. 
.LE. 
.  I.E . 
.LE. 
.LE, 
.  LE . 


RS( 1 ) )  N2=N2+ I 


RS ( 1 ) ) 
RS ( 1 )  ) 
RS( 1 ) ) 
RS( 1) ) 
RS( 1 ) ) 


N3=N3+ 1 
N4  -N4+ 1 
N5=N5+ 1 
N6 -N6  +  1 
N7  =  N7  + 1 


OF  TIMES  THE  TRUE 
CONFIDENCE  INTERVAL*** 
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o  o 


SUBROUTINE  PARES ( N , M , C , THETA , EK , MR , PTH , PER , T ) 

C  INPUT 

C  N=SAMPLE  SIZE  (BEFORE  CENSORING ) ,N= IOO  OR  LESS  AS 
C  DIMENSIONED 

C  SS1=0  IF  SCALE  PARAMETER  THETA  IS  KNOWN 

C  SS1=1  IF  SCALE  PARAMETER  THETA  IS  TO  BE  ESTIMATED 

C  SS2=0  IF  SHAPE  PARAMETER  K  IS  KNOWN 

SS2=1  IF  SHAPE  PARAMETER  K  IS  TO  BE  ESTIMATED 

SS3=0  IF  LOCATION  PARAMETER  C  IS  KNOWN 
C  SS3=1  IF  LOCATION  PARAMETER  C  IS  TO  BE  ESTIMATED 
C  T ( I ) =I-TH  ORDER  STATISTIC  OF  SAMPLE  (I=1,N) 

C  M=NUMBER  OF  OBSERVATIONS  REMAINING  AFTER  CENSORING  N-M 

C  FROM  ABOVE 

C  C ( 1 ) =INITI AL  ESTIMATE  (OR  KNOWN  VALUE)  OF  C 
C  THET A( 1 ) = I N I T I AL  ESTIMATE  (OR  KNOWN  VALUE)  OF  THETA 

C  EK( 1 ) -INITIAL  ESTIMATE  (OR  KNOWN  VALUE)  OF  K 

C  MR=NUMBER  OF  OBSERVATIONS  CENSORED  FROM  BELOW 

C  OUTPUT 

C  N,SS1 ,SS2,SS3,M,C( 1 ) , THETA ( 1 ) ,EK( 1 ) ,MR 
C  — SAME  AS  FOR  INPUT 

C  C ( J ) =E ST I MATE  AFTER  J-I  ITERATIONS 

C  (OR  KNOWN  VALUE)  OF  C 

C  THETA ( J ) =ESTIMATE  AFTER  J-l  ITERATIONS 
C  (OR  KNOWN  VALUE)  OF  THETA 

C  EK( J )=ESTIMATE  AFTER  J-l  ITERATIONS 

C  (OR  KNOWN  VALUE)  OF  K 

C  (MAXIMUM  VALUE  OF  J  AS  PRESENTLY  DIMENSIONED  IS  500) 

C  EL-NATURAL  LOG.  OF  LIKELIHOOD  FOR  C( J ) , THETA ( J ) ,EK( J ) 
DIMENSION  T ( 500 ) , C ( 500 ) , THETA ( 550 ) ,EK( 550 ) ,X(56) ,Y(55) 
SS1=1. 

SS2=1 . 

SS3-0. 

IF ( N ) 66 , 66 , 104 
104  EN=N 

IF (M) 66 , 66 , 1 10 
110  EM-M 
31  ELNM=0 . 

EMR=MR 
MRP=MR+ 1 

33  NM=N-M+1 

DO  34  I =NM , N 
EI  =  I 

34  ELNM  ELNM+ALOG ( El ) 

IF  (MR)  66,35,74 

74  DO  75  1-1 , MR 
EI  =  I 

75  ELNM-  ELNM-AI.OG(EI  ) 

35  DO  30  J=1 ,550 

IF  (J-l)  66,25,37 
37  JJ=J-1 
SK=0 . 

S  L—  0 . 

DO  6  I =MRP,M 

6  SK=SK+ ( T ( I ) — C ( JJ ) ) **EK( J J ) 
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IF(SSl)  7,7,8 

7  THETA ( J ) =THETA( JJ ) 

GO  TO  9 

8  IF  (MR)  66,19,20 

19  THETA ( J ) = ( ( SK+ ( EN-EM ) * ( T(M )-C ( JJ ) **EK( J J ) /EM) 
C**(1./EK(JJ) ) 

GO  TO  9 

20  X( 1 ) =THETA( J J ) 

LS=0 

DO  21  L=1 , 55 
LL=L-1 
LP=L+1 
X( LP ) =X( L ) 

ZRK= ( (T(MRP)-C(JJ) )/X(L) )**EK(JJ) 

Y(L)=-EK( JJ)*(EM-EMR)/X(L)+EK( JJ ) *SK/X( L) **(EK( JJ )  +  l .  ) 
C+EK( J J ) * ( EN-EM )*(T(M)-C( JJ) )**EK( JJ )/X(L)**(EK(JJ)+l. 

C ) -EMR*EK ( J J ) *ZRK*EXP ( -ZRK ) / ( X ( L ) * ( 1 . -EXP ( -ZRK ) ) ) 

IF  (Y(L) )  53,73,54 

53  LS=LS-1 

IF  (LS+L)  58,55,58 

54  LS=LS+1 

IF  (LS-L)  58,56,58 

55  X( LP ) = . 5*X( L ) 

GO  TO  61 

56  X(  LP )  =  1 . 5*X(  L  ) 

GO  TO  61 

58  IF  (Y(L)*Y(LL) )  60,73,59 

59  LL=LL-1 
GO  TO  58 

60  X  (  LP )  =X  (  L  )  +Y  (  L )  *  ( X  ( L )  -X  (  LL )  )  /  ( Y  ( LL )  -Y  (  L ) ) 

61  IF  (ABS(X(LP)-X(L) )-l.E-4)  73,73,21 

21  CONTINUE 

73  THETA (J )=X( LP) 

9  EK( J )-EK( JJ ) 

10  IF  ( SS2  )  12,12,11 

11  DO  17  I=MRP  ,M 

17  SL=SI.+ALOG(T(  I  )-C(  JJ)  ) 

X( 1 )=EK( J ) 

LS=0 

DO  51  I,=  1 ,55 
SLK=0 . 

DO  18  I=MRP ,M 

18  SLK-SLK+ ( ALOG ( T ( 1 )-C ( J J ) ) -ALOG ( THETA ( J ) ) ) * ( T ( I ) -C ( J  J ) ) 
C**X( L ) 

LL=L-1 
LP=L+ 1 
X( LP)-X(L) 

ZRK-  (  ( T  ( MRP )  -C  (  J  J  )  )  /THETA  ( J  )  )**X(  I.) 

Y ( L )  —  { EM-EMR ) * ( 1 ./XI I.) -ALOG ( THETA ( J ) I ) +SL-SLK/THETA( J ) 
C**X  ( I. )  +  ( EN-EM )  *  I  ALOG  I  THETA  (  J  )  )  -ALOG  (  T  ( M  )  -C  ( J  J  )  )  )  *  (  T I  M  ) 
C-C( JJ) )**X( L ) /THETA ( J ) **X(  I, ) +EMH*ZKK*  <  ALOG ( ZRK ) /X( L ) ) 
C*EXP(-ZRK)/( 1 .-EXP(-ZRK) ) 

IF ( Y( L ) )  43,52,44 
43  LS-LS-i 
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IFUS  +  L)  47,45,47 

44  LS=LS+1 

IF  (LS-L)  47,46,47  - 

45  X(LP)=.5*X(L) 

GO  TO  50  . 

46  X(  I;P)  =  1 . 5*X(  L ) 

GO  TO  50 

47  IF  (Y(L)*Y(LL) )  49,52,48 

48  LL=LL-1 
GO  TO  47 

49  X( LP)=X(L) +Y( L) *(X( L)-X( LL ) ) /( Y( LL)-Y( L) ) 

50  IF  ( ABS (X( LP)-X(L) )-l .E-4)  52,52,51 

51  EK( J )=X( LP) 

12  C(J)=C(JJ) 

62  IF  ( SS3 )  25,25,14 

14  IF  (l.-EK(J))  16,78,78 

78  IF  (SS1+SS2)  57,57,16 
16  X(1)=C(J) 

LS=0 

DO  23  L=1 ,55 
SK1=0. 

SR=0 . 

DO  15  I =MRP , M 

SK1=SK1+(T(I ) —X ( L ) )**(EK( J )-l . ) 

15  SR=SR+ 1 ./(T(I)-X(LL) ) 

LL=L-1 

LP=L+1 

X(LP)=X(L) 

ZRK=( ( T ( MRP )-X( L ) ) /THETA (J ) ) **EK( J ) 

Y ( L ) = ( 1 . -EK( J ) ) *SR+EK( J ) * ( SKI + ( EN-EM ) * ( T ( M ) -X ( L ) ) 

C** ( EK( J ) -1 . ) )/THETA(J)**EK(J ) -EMR*EK( J ) *ZHK( EXP( -2RK ) 
C/( (T(MRP)-X(L) )*( 1 .-EXP(-ZHK) ) ) 

IF  (Y(L) )  39,24,40 

39  LS=LS-1 

IF  (LS+L)  70,41,70 

40  LS=LS+ 1 

IF  (LS-L)  70,42,70 

41  X(  LP )  =  . 5*X( L ) 

GO  TO  22 

42  X(LP)-.5*X(L)+.5*T( 1 ) 

GO  TO  22 

70  IF  (Y(L)*Y(LI.)  )  72,24,71 

71  LL=LL-1 
GO  TO  70 

72  X(LP)=X(L)+Y(L)*(X(L)-X(LL) )/(Y( LL)-Y(L) ) 

22  IF  (ABS(X(l.P)-X(L)  )-l.E-4)  24,24,23 

23  CONTINUE 

24  C ( J ) =X( LP ) 

GO  TO  25 

57  C(J)=T( 1 ) 

25  IF  (MR)  66,38,69 
38  DO  63  I - 1 , M 

IF  (C( J)+l .E-4-T( I ) )  68,67,67 
67  MR-MR+ 1 
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63  C ( 1 >  =T ( 1 ) 

68  IF (MR)  66,69,31 

69  SK=0. 

SL=0. 

DO  36  I =MRP ,M 
SK=SK+( T (I)— C(J))**EK(J) 

36  SL=SL+ALOG ( T ( 1 ) -C ( J ) ) 

ZRK= ( (T(MRP)-C ( J ) ) /THETA (J ) ) **EK( J ) 

EL=ELNM+ ( EM-EMR ) * ( ALOG ( EK  ( J  ) )-EK( J ) *ALOG ( THETA ( J ) ) )+ 
C (EK( J )-l . ) *SL-( SK+ (EN-EM) * ( T(M )-C ( J ) ) **EK( J ))/( THETA 
C( J )**EK( J ) )+EMR*ALOG( 1 .-EXP(-ZRK) ) 

150  IF( J-3)  30,27,27 

27  IF  ( ABS (C(J)-C(JJ) ) — 1 .E — 4)  28,28,30 

28  IF  ( ABS (THETA (J )— THETA( JJ ) )-l .E-4)  29,29,30 

29  IF(ABS(EK(J)-EK( J J ) )- l .E-4) 126 , 126 ,30 

30  CONTINUE 
126  PTH=THETA( J ) 

PEK=EK( J ) 

GO  TO  140 
66  PRINT  135 

135  FORMAT ( 1H  ,20HAL  SAMPLES  CENSORED,/) 

PEK=0 . 

PTH-0 . 

140  CONTINUE 
RETURN 
END 


59 


VITA 


James  Ward  Johnston,  Jr.,  was  born  on  28  September 
1950  in  San  Diego,  California.  He  graduated  from  Point 
Loma  High  School  in  San  Diego  in  1968  and  attended  San  Diego 
State  University,  where  upon  graduation  in  1973  he  received 
a  bachelor 's  degree  in  Mathematics.  He  then  worked  as 
a  cadastrial  draftsman  for  the  County  of  San  Diego  for  a 
period  of  one  year.  He  then  completed  training  at  Officers 
Training  School  at  Lackland  AFB  TX,  and  was  commissioned  a 
second  lieutenant  in  the  United  States  Air  Force.  His 
initial  assignment  was  to  the  E-4  System  Program  Office  at 
Hanscom  AFB  MA,  where  he  served  as  a  Computer  Design  Engi¬ 
neer.  After  approximately  four  years  at  Hanscom  AFB,  he 
entered  the  School  of  Engineering,  Air  Force  Institute  of 
Technology,  in  June  1979. 


Permanent  address:  3235  Garrison  Street 

San  Diego,  California  92106 


60 


_ UNCLASSIFIED _ _ 

SECURITY  CLASSIFICATION  OF  THIS  PACE  (When  Halt.  hiitunJ) 

REPORT  DOCUMENTATION  PAGE 

1  REPORT  NUMBER  |?  GOVT  ACCESSION  NO 

AFIT/GOR/OS/8QD  -5 _ /tPV/9& 

4.  TITLE  (end  Subtitle)  * 

A  MODIFIED  DOUBLE  MONTE  CARLO  TECHNIQUE  TO 
APPROXIMATE  RELIABILITY  CONFIDENCE  LIMITS 
OF  SYSTEM  WITH  COMPONENTS  CHARACTERIZED  BY 
THE  WE I BULL  DISTRIBUTION 

7  authorc*; 


KK  Al>  INSTRUCT  IONS 
HI  I- OKI-.  COMITT  I  INI.  I- 1 1  KM 
recipient's  catm.ii‘,  niimher 


TYPE  Of  REPORT  4  ptRIOD  COYEHED 

MS  Tiles  i_s 

PERFORMING  OTG.  REPORT  NUMBER 
CON  T  R  AC  T  OfT G~R  A  N  T  N  U  M  BE  R<  v> 


Janies  W.  Johnston,  Jr.,  Captain,  USAF 


9  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS  0 

Air  Force  Institute  of  Technology  (AEIT/ENl 
Wright-Patterson  AFB  OH  45433 


o  CONTROLLi  iG  OFFICE  NAHIl  AND  ADDRESS  Z 

Air  Force  Institute  of  Technology  (AFIT/ENi-,- 


PROGRAM  EL  EMEN  T  PROJECT  TASK 
AREA  &  WORK  UNIT  NUMBERS 


REPORT  DATE 

December  1 980 

TTuMBFR  OF  PAG E S 


U  MONITORING"  AGENCY  NAME  4  AOORESSfl/  ililforeiil  Irom  Cunlrolllnn  Oltit  a) 


_ _ 60  _ 

is  SECURITY  CLASS  (ol  ihla  r»p,>rll 


UNCLASSIFIED 

'•.«  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 


16  DISTRIBUTION  STATEMENT  -of  thla  Report) 


Approved  for  public  release;  distribution  unlimited 


SUPPLEMENTARY  NOTES 


Approved  for  Public  Release,  AFR  190-17 


LAUREL  A.  LAMP  El. A ,  2d  Lt  ,  USA!-' 
Upputy  I)  i  i  1‘iTor  ,  l’ul>lic  Allairs 

19  KEY  WORDS  (Continue  on  se  aide  il  ne«  mut  Identity  by  black  number) 

Weibull  Distribution  Reliability 
Double  Monte  Carle 
System  Reliability 
Confidence  Limits 


06  JAN  1981 


20  ABSTRACT  (Continue  on  r*vci  ia  ii.lt  (/  nn  n-ioarv 

~A  modified  Double  Monte 
mi  nation  of  lower  confident 
component  test  data.  It  is 
consists  of  failure  times 
known  two-parame ter  Weibull 
failure  times  are  randomly 
scale  parameters  of  the  dis 


end  tdentiU  fi>  blink  number) 

Carlo  technique  is  developed  for  deter- 
e  limits  of  system  reliability  based  on 
assumed  that  the  components  test  data 
which  are  distributed  according  to  a 
probability  distribution.  These 

generated  using  the  true  shape  and  _ 

tribution.  Max  I  mum- I i ke I i hood 


00 


FORM 

i  l  AN  ti 


147  3 


I  M  »  V  *  S  IS 


UNC LASS  1  FI KD 

|UiTY  "  l  A^ilFti  AT  (ON  f  THIS  »  'Gf  .  • 


i  ereil 


L  i 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGEC*?!./!  Dml.  I  nt.rtd) 

- ; — — - - 

| estimators  are  twice  obtained  lor  the  shape  and  scale  parameters 

and  then  substituted  into  the  reliability  equation  to  obtain  an 
-stimator  for  the  component  reliability.  A  s' i  veil  number  of  i  hose 
estimators  are »ob ta i ned  and  n  ■  d  to  form  an  empirical  distribution 
of  reliabilities  for  each  component.  A  given  niimher  ol  samples 
from  this  distribution  are  used  to  calculate  various  system 
reliabilities.  Since  the  trie-  system  reliability  is  known,  it  can 
be  determined  if  a  given  confidence  interval  contain  the  i  i-ii< 
number,  hence  giving  you  a  method  of  validating  any  confidence 
interva 1 . 

r 

\ 


